# ******** Study the Simulation Results From the Model *******
### Note: the default is that when you call this file, the path is the current file path. "FigurePath" and "TablePath" are paths for figure and table outputs. 

# ------------  Pre-Process the Model Simulated Data -----------------
# Then call the module that processes
option="main_paper"   #pre-process for the main paper.  Otherwise, select "general_belief", and then the file will process results for general beliefs, and those will be in the appendix.
crisis_def = "bank_equity_change"  # Define crisis as based on bank equity changing below a threshold  
source("Codes_Model/preprocessing_simulation_data.R")
# Then select whether to generate tables to files. 
gen_tables = TRUE

#==================== Table 4 Panel A: Main regression results ===================
print("************ Generating Table 4 Panel A ************ ")
reg_output_and_credit_spread <- function( TB ){
  TB[ , crisis_interact_credit_spread_change :=  credit_spread_change*crisis]
  TB[ , crisis_interact_credit_spread_control :=  credit_spread_last_year*crisis]
  TB[  , crisis_interact_bank_credit:= bank_credit_GDP*crisis ]
  regs=list()
  regs[[1]] = lm( 100*GDP_3Y_change ~  crisis  , data=TB )
  regs[[2]] = lm( 100*GDP_3Y_change ~  crisis_interact_bank_credit  , data=TB )
  return(regs)
}
latex_reg_output_and_credit_spread_all <- function(regs, regs_data, dependent_var_name="GDP Growth"){
  N_columns = length(regs)+length(regs_data)
  results=capture.output(stargazer( c(regs, regs_data) , df = FALSE, digits = 2, omit.stat = c("rsq", "f", "ser","n","adj.rsq") , omit=c("Constant"), star.cutoffs = c(0, 0, 0),
                                    se = list( 0, 0, 0, 0, 0)  ,
                                    title=paste0("Crisis Severity, Credit Spread and Bank Credit (using equity crash to define a crisis)"),
                                    label = paste0("Severity regression on equity crashes")  ,
                                    dep.var.caption= paste0("\\textit{Dependent variable: ", dependent_var_name, " from $t$ to $t+3$}") ,
                                    dep.var.labels =  rep("",N_columns) , 
                                    column.labels   = c("Bayesian","Diagnostic", "Data"),
                                    covariate.labels = c("crisis$_t$", 
                                                         "$(\\frac{\\text{bank credit}}{\\text{GDP}})_t*$crisis$_t$" ),
                                    column.separate = c(2,2,1), 
                                    add.lines =  list(  c( "Observations", rep("",length(regs)), 2548  )  ),
                                    notes="To be replaced", notes.append=FALSE, notes.align="l",
                                    column.sep.width = "16pt", no.space = TRUE
  ))
  results[12] = "\\\\[-1.8ex]"
  results[length(results)-2] =  paste0("\\multicolumn{", N_columns+1 ,"}{p{0.9\\textwidth}}{\\footnotesize 
                                       \\textit{Note}: Model and data regressions are normalized so that the coefficients reflect the impact of one sigma change in spreads, and bank credit/GDP. 
                                       The coefficient in column (5) is from Table I (column 4) of \\cite{baron2021banking}. 
                                       }")
  writeLines(results, paste0(TablePath,"crisis_aftermath/",gsub(" ", "_", dependent_var_name, fixed = TRUE), "_all_equity_defined_crisis.tex" ) )
  print(results)
}
regs_benchmark = reg_output_and_credit_spread(TB_benchmark)
regs_Bayesian = reg_output_and_credit_spread(TB_Bayesian)
regs_Diagnostic= reg_output_and_credit_spread(TB_Diagnostic)
regs_data = regs_benchmark[[1]]
regs_data$coefficients[2]=-4.5
if(gen_tables){
  regs = c(regs_Bayesian, regs_Diagnostic )
  latex_reg_output_and_credit_spread_all( regs, list(regs_data) )
}

#  ================== Table 4 Panel B: Predicting Equity Crashes  ==================
print("************ Generating Table 4 Panel B ************ ")
reg_crises_predictability <- function(TB, auxiliary=FALSE){
  TB_yearly = TB
  cutoff = 0.5
  TB_yearly[ , froth := credit_spread<=quantile(credit_spread,cutoff) ]
  TB_yearly[ , crisis_next_year :=  c(  crisis[2:nrow(TB_yearly)]>0,  NA  ) ]
  TB_yearly[ , crisis_next_5_year :=  c(  rollsum(crisis,5)>0, rep(0,4) ) ]
  TB_yearly[ , high_credit := bank_credit_GDP>=quantile(bank_credit_GDP,1-cutoff)  ]
  TB_yearly[ , high_credit_past_5y := c( rep(0,4), rollmean(high_credit,5 ) )>0 ]
  TB_yearly[ , capital_crash := capital_excess_return<quantile(capital_excess_return,0.05) ]
  marginal_effect <- function(reg){
    summary_results = summary(margins(reg))
    return(as.numeric(summary_results[2])*100)
  }
  
  regs = list()
  regs[[1]] = glm( crisis_next_5_year ~ bank_credit_GDP, family = binomial(link = "probit"), data = TB_yearly)
  regs[[1]]$coefficients[2] = marginal_effect(regs[[1]])
  return(regs)
}
regs_data = reg_crises_predictability(TB_Bayesian)
regs_data[[1]]$coefficients[2] = 5.4

latex_reg_crises_predictability_all <- function(regs){
  regs_output = c(regs, regs_data)
  results=capture.output(stargazer( regs_output , df = FALSE, digits = 2, omit.stat = c("rsq", "f", "ser","n","adj.rsq","LL","aic") , omit=c("Constant"), star.cutoffs = c(0, 0, 0),
                                    se = list( 0, 0, 0)  ,
                                    title=paste0("Predicting Equity Crashes"),
                                    label = paste0("prediction of equity crashes")  ,
                                    dep.var.caption= paste0("\\textit{Dependent variable: equity crash$_{t+1\\text{ to }t+5}$}") ,
                                    dep.var.labels =  rep("",3) , 
                                    column.labels   = c("Bayesian", "Diagnostic", "Data"),
                                    covariate.labels = c("$(\\frac{\\text{bank credit}}{\\text{GDP}})_t$"),
                                    column.separate = c(1,1,1), 
                                    add.lines =  list(  c( "Observations", rep("",2), c(316)  )  ),
                                    notes="To be replaced", notes.append=FALSE, notes.align="l",
                                    column.sep.width = "16pt", 
                                    no.space = TRUE
  ))
  results[12] = "\\\\[-1.8ex]"
  results[length(results)-2] =  "\\multicolumn{4}{p{0.75\\textwidth}}{\\footnotesize \\textit{Note}: 
  The coefficient on Bank Credit/GDP is the sensitivity of crisis probability (\\%) to a one standard deviation increase in bank credit/GDP. 
  The data regression is from Table III (column 7) of \\cite{baron2017credit}.
}  "
  writeLines(results, paste0(TablePath,"predictability/predicting_crisis_all_equity_defined_crisis.tex" ) )
  print(results)
}

regs_benchmark = reg_crises_predictability(TB_benchmark)
regs_Bayesian = reg_crises_predictability(TB_Bayesian)
regs_Diagnostic = reg_crises_predictability(TB_Diagnostic)
regs = c( regs_benchmark, regs_Bayesian, regs_Diagnostic )
stargazer( regs, type="text", df = FALSE, digits = 2, omit.stat = c("rsq", "f", "ser") , star.cutoffs = c(0, 0, 0 ), 
           column.labels   = c("Static Belief", "Bayesian", "Diagnostic"), column.separate = rep(length(regs)/3,3), 
           dep.var.labels.include = FALSE, 
           covariate.labels = c("bank credit/GDP"),  dep.var.caption = "Dep var: equity crash in the next 5 years" )

if(gen_tables){
  latex_reg_crises_predictability_all( c(regs_Bayesian,regs_Diagnostic) )
}

